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We have observed atmospheric gamma rays from 30 GeV to 8 TeV, using emulsion cham- 
bers at balloon altitudes, accumulating the largest total exposure in this energy range to date, 
SQ.T ~6.66 m 2 sr day. At very high altitudes, with residual overburden only a few g cm -2 , at- 
mospheric gamma rays are mainly produced by a single interaction of primary cosmic rays with 
overlying atmospheric nuclei. Thus, we can use these gamma rays to study the spectrum of primary 
cosmic rays and their products in the atmosphere. From the observed atmospheric gamma ray 
spectrum, we deconvolved the primary cosmic-ray proton spectrum, assuming appropriate hadronic 
interaction models. Our deconvolved proton spectrum covers the energy range from 200 GeV to 
50 TeV, which fills a gap in the currently available primary cosmic-ray proton spectra. We also 
estimated the atmospheric muon spectrum above 30 GeV at high altitude from our gamma-ray 
spectrum, almost without reference to the primary cosmic rays, and compared the estimated flux 
with direct muon observations below 10 GeV. 

PACS numbers: 92.60.hx, 13.85.Tp 



I. INTRODUCTION 

The primary cosmic-ray proton spectrum is one of the 
most important quantities needed to interpret cosmic-ray 
phenomena inside the atmosphere. Many observations of 
primary protons have been performed since the discovery 
of cosmic rays. Recently, proton flux measurements have 
attracted renewed attention, since they are needed to 
precisely estimate the absolute flux of atmospheric neu- 
trinos, providing critical input for analyses of neutrino 
oscillations being performed with Super-Kamiokande 
and similar detectors. In particular, atmospheric neutri- 
nos with higher energies, ~ 100 GeV, which are detected 
as upward through-going muons, are mainly produced 
from the interaction of ~ 1 TeV primary protons in the 
atmosphere. Similarly important are observations of the 
absolute flux of muons, which are the partners of neutri- 
nos, and thus provide a check on atmospheric neutrino 
flux estimates. Since atmospheric gamma rays at high 
altitude, with residual overburden of a few g cm -2 , are 
mainly produced by a single interaction of primary cos- 
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mic rays with overlying atmospheric nuclei, we can use 
these gamma rays to study the spectrum of primary cos- 
mic rays and their products in the atmosphere. 

We have been observing primary cosmic ray electrons 
at balloon altitudes using emulsion chambers for many 
years, and successfully obtained the energy spectrum of 
electrons in the energy range from 30GeV to 3TeV 0,13] 
In the course of these electron observations, we have si- 
multaneously observed atmospheric gamma rays. In par- 
ticular, since gamma rays produce electromagnetic show- 
ers just like electrons, emulsion chambers are not only 
suitable to observe the gamma rays in the atmosphere, 
but measurement of the gamma ray showers is a neces- 
sary part of the electron analysis, required to check and 
calibrate the performance of the emulsion chambers. 

Atmospheric gamma rays have been observed by the 
BETS group in the energy range from a few GeV to sev- 
eral 10 GeV at mountain altitude (2.77 km) and at bal- 
loon altitudes (15 — 32 km) They found that hadronic 
interaction simulation codes such as the Lund group's 
Fritiof V7.02, or Dpmjet3 give results fairly consistent 
with the observed energy spectra and altitude variation 
of the gamma-ray flux, within the accuracy of their mea- 
surements. The JACEE group also observed high-energy 
gamma rays in the 3 — 30 TeV range at overburden 
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-5.5 g cm" 2 for 120 - 150 hr [|. 

Atmospheric neutrinos, muons, and gamma rays are 
decay products of pions produced in the interaction of 
primary cosmic rays with atmospheric nuclei: 

p + N -» TT^/^ + Z^^) (1) 
\e ± + v e (v e ) + v^ivy) 

-> 7T° -> 2 7 

In addition to pions, if and 77 mesons are produced 
in the primary interactions, and also decay into muons 
and gamma rays, making minor contributions to the net 
fluxes. 

Atmospheric gamma rays are mainly the decay prod- 
ucts of 7T°. When we observe atmospheric gamma rays 
at very high altitudes, corresponding to residual over- 
burdens of several g cm -2 , the thickness of atmosphere 
above the detector is —0.05 nuclear mean free path or 
—0.1 radiation lengths (r.L). Thus, hadronic cascades 
and electromagnetic showers cannot develop significantly 
in the atmosphere before reaching the detector. The ob- 
served gamma rays are almost all produced by the first 
interaction of primary cosmic rays (protons, helium and 
heavier nuclei) with atmospheric nuclei. Therefore, us- 
ing the appropriate hadronic interaction models, we can 
reliably deconvolve the primary proton flux from the ob- 
served spectrum of atmospheric gamma rays produced 
by 7T°. This can be done in a semi- analytic way, because 
of the approximate scaling nature of the pion production 
cross section. In this deconvolution, we include the ef- 
fects of contributions from primary helium and heavier 
nuclei, and also the minor contribution of gamma rays 
from decays of rj and K mesons. 

The proton spectra have been measured with two 
kinds of detectors: magnetic spectrometers and calorime- 
ters. Although magnetic spectrometers have excellent 
energy resolution, their maximum observable energy is 
limited by the maximum detectable rigidity (MDR) to 
-1 TeV. The BESS, AMS-01 and CAPRICE instruments 
are typical detectors using magnetic spectrometers, and 
have precisely measured primary proton spectra up to 
- 100 GeV 0, Sll- BESS and AMS-01 give results con- 
sistent with each other within an error of a few %, while 
CAPRICE reports a 10 % lower proton flux. Recently, 
the BESS-TeV spectrometer, with improved MDR [<|, 
observed primary proton spectra up to 540 GeV. As for 
calorimeters, electronic detectors using scintillators, as 
well as passive emulsion chambers have been used to ob- 
serve the proton spectrum in the higher energy region. 
The first measurements of the proton spectrum in the 
TeV region were made by the PROTON satellite-borne 
calorimeters |10| . These observations were followed by an 
ionization spectrometer flown at balloon altitudes whose 
energy range extend from 50 GeV to 2 TeV 11]. The 
proton spectrum in the energy range from 10 TeV to 
1000 TeV has been measured by the JACEE Collabora- 
tion and RUNJOB Collaborations using emulsion cham- 



ber detectors [2 E3- 13 and references therein]. Re- 
cently, ATIC, which is a calorimeter with BGO scintil- 
lators, observed primary proton spectrum in 30 GeV- 
50 TeV ^3- However, more accurate measurements of 
the proton spectrum are still required in the energy range 
from 100 GeV to 10 TeV. 

In detectors that directly observe primary cosmic rays, 
such as JACEE and RUNJOB, the energies of primary 
protons are estimated from electromagnetic showers pro- 
duced by nuclear interactions within the emulsion cham- 
bers. Since the electromagnetic shower developed in 
the chamber is actually initiated by multiple secondary 
gamma rays, and its structure is affected by their emis- 
sion angles, energy estimation for the primary particle is 
complicated. On the other hand, in the present work, we 
derive the energy of primary protons from observed elec- 
tromagnetic showers initiated by one single gamma ray, 
so that the estimation of proton fluxes from atmospheric 
gamma rays is simpler and therefore more reliable. 

As a first approximation, the number of charged pions 
produced in hadronic interactions is almost two times of 
neutral pions. Hence we can estimate the production rate 
of charged pions from atmospheric gamma rays. Since 
the muons are mostly produced from the decay of charged 
pions, we can also directly estimate the muon flux with- 
out reference to the primary cosmic ray flux, correcting 
for the minor contribution of rj and K mesons. More pre- 
cisely, widely accepted hadronic interaction models are 
used to derive the muon spectrum from the gamma-ray 
spectrum. 

The MASS, CAPRICE, HEAT and BESS group have 
performed direct atmospheric muon observations at var- 
ious balloon altitudes [H E3, El [H H|. The muon 
spectrum at various atmospheric depths gives important 
information to check the reliability of calculations per- 
formed to estimate the neutrino flux [21|. However, the 
muon flux at large depths depends on the detailed struc- 
ture of atmospheric density as a function of altitude, 
which is subject to seasonal variation and meteorologi- 
cal conditions. This brings complications and ambigui- 
ties into the muon flux estimation. On the other hand, 
the flux of atmospheric gamma rays at stratospheric alti- 
tudes depends only on the residual atmospheric overbur- 
den, which can be estimated much less ambiguously than 
the muon flux. Hence, an accurate measurement of the 
gamma-ray flux can be also used to calibrate the neu- 
trino flux calculations, because of the relation between 
gamma-ray and neutrino production in the atmosphere 
as shown in the process (Q. 

In this paper, we present measurements of the at- 
mospheric gamma-ray spectrum due to hadronic inter- 
actions in the energy range of 30 GeV to 8 TeV, ob- 
served at balloon altitudes with emulsion chambers. We 
also present results from the deconvolution of primary 
cosmic-ray proton and high altitude muon spectra from 
our gamma-ray observations. 
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II. ATMOSPHERIC GAMMA-RAY 
OBSERVATION 

A. Balloon observations 

We have observed primary electrons with balloon- 
borne emulsion chambers in many flights between 1968 
and 2001 @, 0, l^- Simultaneously, we have also ob- 
served atmospheric gamma-rays to check the perfor- 
mance of the emulsion chambers in each balloon experi- 
ment. The pressure altitude records for each flight corre- 
spond to residual atmospheric overburdens in the range 
from 4.0 g cm~ 2 to 9.4 g cm~ 2 . The total cumulative 
effective exposure SflT for gamma rays is 6.66m 2 -sr-day, 
which is larger than any other atmospheric gamma-ray 
observations performed at balloon altitudes in the energy 
range of 30 GeV to 8 TeV. In Table [I] we summarize the 
series of experiments since 1968. Nishimura et al. (1980) 
have reported the spectra of primary electrons and at- 
mospheric gamma rays for the 1968 — 1976 observations 
in our previous work 

B. Detector 

Emulsion chambers consist of nuclear emulsion plates 
and lead plates, which are stacked alternately. Nuclear 
emulsion plates sample the development of electromag- 
netic showers produced in the lead plates. X-ray films are 
also inserted to allow rapid, naked-eye scanning for high 
energy showers, which produce dark spots in the film. 
The threshold energy for shower detection depends on the 
background and the sensitivity of the films; some details 
of the detection threshold for various films are described 
later. Figure ^ shows a typical emulsion chamber con- 
figuration. Detailed configurations and performance are 
described in Nishimura et al. (1980) 2] and Kobayashi 
et al. (1999) 0. 

In emulsion chambers, it is possible to measure the 
location of shower tracks in each emulsion plate with a 
precision of ~1 fim. Because of this high position resolu- 
tion, we can inspect the shower starting points in detail 
and unambiguously distinguish showers due to electrons, 
gamma rays, and other hadronic interaction events 0- 
By inspecting various specific features of those events, 
the rejection power for protons misidentified as electron 
candidates is found to be as large as 10 5 Q. 

We measure the shower particles within a circle of 
100 fiui radius from shower axis, which means that we se- 
lect the shower particles with higher energy, which have 
suffered less multiple scattering in the chamber. Hence, 
the number of the shower particles drops off faster than 
for all shower particles, and the shower maximum ap- 
pears in ~6 r.l. for 1 TeV electrons, while the maximum 
of the total number of shower particles appears in ^12 
r.l. for 1 TeV electrons. This means that we can de- 
termine the energy of higher energy incident electrons 
with thinner detector. The typical size and thickness of 



the detector are 40 cm x 50 cm, and 8 cm (~9 r.l.), 
respectively. Thus the emulsion chamber has the advan- 
tage of a large effective area combined with a wide field 
of view compared with other detectors. Because of the 
simple configuration of the detector, we can estimate the 
geometrical factor (SQ) very accurately, a task which is 
difficult for some electronic detectors. 

For electron observations, the effective geometrical fac- 
tor is given by 

SQ e = 2irSr] [ cos6sin0d9 = 7rS?7sin 2 6>o, (2) 
Jo 

where 0q is the upper limit of incident angles and rj is 
the efficiency of events that pass through the top and 
bottom emulsion plates, the so-called "edge effect". On 
the other hand, for atmospheric gamma-ray observations, 
the effective geometrical factor is given by 

Sn 7 = 2ttSt] / sin<9d<9 = 27rS*77(l-cos<9o). (3) 
Jo 

Here, we corrected 5f2 7 to obtain the vertical flux, be- 
cause the gamma rays have a l/cos# enhancement rel- 
ative to isotropic primary cosmic rays, proportional to 
the path length in the overlying atmosphere at balloon 
altitudes. Figure |3 shows the zenith angle distribution 
observed with emulsion chambers, compared with expec- 
tation calculated taking into account small corrections 
due to the elongation of the path length in the overly- 
ing air for large zenith angles, which affect the amount 
of absorption of the electrons and gamma rays. In the 
typical case of 6 = 60° and S = 0.40x0.50 m 2 , Sft 7 is 
0.52 m 2 sr with rj of 0.82 for the chamber thickness of 
8.0 cm. For the comparison of emulsion chambers with 
other detectors, we summarize the effective geometrical 
factor SQ e with efficiency for primary electron experi- 
ments in Table [H] We can see how the emulsion cham- 
ber is efficient for observation of low-flux components of 
cosmic rays such as primary electrons and atmospheric 
gamma rays. 

1. Scanning method 

Since high energy electro-magnetic showers above sev- 
eral 100 GeV leave dark spots on X-ray films, we can 
detect these showers with the naked eye by scanning the 
X-ray films. We locate the corresponding tracks in the 
adjacent emulsion plate using microscopes, trace them 
back through the stack to the cascade starting point, 
and identify the incoming particle. The validity of par- 
ticle identification in our experiments is also checked by 
comparison with the expected zenith angle distributions 
and shower starting point distributions, as described in 
our previous work |2j,|3j. We picked up events with zenith 
angle less than 60°. The detection threshold of the X-ray 
film, although it depends on the accumulated background 
tracks and fog in the film, is ^500 GeV for Sakura type-N 
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TABLE I: List of balloon flights 
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* Effective ST2 7 T for the gamma-ray observations. 




FIG. 1: Typical configuration of the emulsion chamber in 
cross-sectional drawing from side view. 



X-ray film used before 1984, -800 GeV for Fuji #200 X- 
ray film and —200 GeV for screen type X-ray films such as 
Fuji G8-RXO, G12-RXO, HR8-HA30, HR12-HA30 used 
on and after 1984 

To detect the electro-magnetic showers below a few 
hundred GeV, the emulsion plates had to be directly 
scanned using microscopes. Although the emulsion 
chamber is the only detector which has succeeded in ob- 
serving cosmic-ray electrons and gamma rays in the TeV 
region, it was difficult to observe electrons or gamma-rays 
below a few hundred GeV because of the tedious work 



FIG. 2: Zenith angle distribution for atmospheric gamma 
rays beyond 400 GeV. 



required in manual microscope scanning of the emulsion 
plates. In order to overcome this difficulty tracing low en- 
ergy showers, in 2001 we started to apply the automatic 
scanning method which has been successfully developed 
by the Nagoya group to observe v T events with nuclear 
emulsions 29]. A fully automatic emulsion scanning 
system developed by the Nagoya group, called "Track 
Selector", consists of a motor-driven microscope XYZ 
stage with a controller, a fast read-out CCD and special 
hardware for 3D image processing,. It takes the images 
of 16 "slices" , or layers within an emulsion plate, each 
of 512x512 pixels, corresponding to a field of view of 
147x106 fim 2 , and then outputs position and angle data 
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TABLE II: Examples of effective geometrical factors with 
efficiency of the electron detectors. 



Detector 


Sfl e (m*sr) 


reference 


T.R.D. 


0.11 


|23J 


MASS-91 


0.018 


f24] 


CAPRICE94 


0.016 
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HEAT 


0.012 


f26] 


BETS 


0.032 


[27] 


AMS-01 


0.10 


\7] 


ECC 


0.38* 


® 



ST2 e in the typical emulsion chamber. 



fluctuations of the number of shower particles, while the 
last term represents fluctuations due to "punch-through" , 
shower particles escaping from the finite thickness of the 
detector. The adjustable coefficients of a, b, and c are 
derived for each chamber using the results of the Monte 
Carlo simulations. Figure El shows some examples of 
the energy dependence of energy resolution for incident 
gamma rays from the simulations, whose showers start 
from pair electrons. The fitted functions of Q are also 
plotted. The coefficients of a, b, and c typically have the 
values of a~\2 - 13 %, 6-5 - 10 %, and c~2 - 3 %. 



for each track found in the emulsion plate. It takes ap- 
proximately 16 hr to perform general scanning over an 
area of 1 cm 2 . We scan the emulsion plate at 3 r.l. depth 
with the Track Selector, and read out all tracks with 
tan6> < 0.3. Shower candidates in the track data are 
identified with off-line programs. Selected shower can- 
didates are traced back to the top plate of the chamber 
by searching for tracks having the predicted angle and 
position. We carefully inspect the starting behavior of 
a shower to see whether it was initiated by an electron, 
gamma ray or hadron. For the balloon detectors in 2001, 
we successfully observed atmospheric gamma rays down 
to 30GeV using this automatic scanning system. 



2. Energy determination 
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Shower energy was determined by comparing the num- 
ber of shower tracks at various depths with the theo- 
retical transition curves, and fitting to the integrated 
track length, used to estimate total ionization in the 
cascade. As the chamber structure is slightly different 
for each flight, we calculated the shower development 
for each chamber using a Monte Carlo simulation code 
called EPICS H3- Results calculated using the EPICS 
code were confirmed by emulsion chambers exposed to 
FNAL electron beams in which showers of 100 GeV 
electrons were re-analyzed. Figure |3| 0] show longitu- 
dinal development of the average number of shower elec- 
trons, and energy distribution of the FNAL experimental 
and simulated data. The simulations well represent the 
experimental data, and the determined energy with the 
simulation for 100 GeV electrons is consistent with the 
experiment within a precision of a few %. The energy 
resolution is 12 % at 100 GeV, as shown in Fig.0] The 
energy resolutions for each emulsion chamber are well 
represented by the form of 



JL = [ 2/ ^0 

E [a MOOGeV 



)- 1 +6 2 + c 2 ( 



E 



lOOGeV 



1/2 



(4) 



where E is the incident gamma-ray energy and a is the 
standard deviation of energy determination. The first 
term in right-hand side root represents statistics-related 



FIG. 3: Longitudinal development of the averaged number of 
shower electrons within a radius of 100 /im. Error bars show 
the uncertainties of the means. Theoretical curve is taken 
from Nishimura et al. (1980) 0. 



C. Atmospheric gamma-ray spectrum 

We observed atmospheric gamma rays at each balloon 
altitude, and derived the vertical gamma-ray spectrum 
normalized at 4.0 g cm~ 2 originated from hadronic inter- 
actions using the following formula: 



ME) = ( 



v. 



Sfl~fT AEC c ffC cn h 



— Cbrem)'Calt 



(m-W^GeV" 1 ). (5) 



Here iV 7 is the number of gamma-ray events, C c ff is 
gamma-ray detection efficiency, C cn h is enhancement fac- 
tor due to the energy resolution as described below, 
Cbrem is bremsstrahlung gamma-ray flux from primary 
electrons, and C a i t is an altitude conversion factor to 
4.0 g cm" 2 . STLy is the geometrical factor to obtain the 
vertical flux. 
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FIG. 4: Energy distributions of the simulated and experi- 
mental data for 100 GeV electron beams at FNAL. 
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FIG. 5: Examples of the energy dependence of energy reso- 
lutions with the emulsion chambers for gamma rays from the 
simulations. The dash lines show the fitted function of P). 



Gamma-ray detection efficiency C c ff in the chamber is 
given by 



Ceff = 1 - exp(-(T r c ), 



(6) 



where <7n = 0.7733 is the probability of pair creation in 
one radiation length, T c is the threshold depth of shower 
starting points for gamma rays. C e s ranges from 0.902 to 
0.955 for our observations, corresponding to T c = 3.0 — 
4.0 r.l. in the threshold depth. 

The uncertainty of the energy determination has the 
effect of enhancing the absolute flux of gamma rays, in 
particular, for the steep power-law spectrum. We de- 
rived the enhancement factor C en h due to the energy res- 
olution for each chamber, which ranges from 1.00 to 1.06 



depending on gamma-ray energies and emulsion chamber 
structures (see appendix lAl in detail). 

In this atmospheric gamma-ray spectrum, we sub- 
tracted the bremsstrahlung gamma rays produced by pri- 
mary cosmic-ray electrons. The observed primary elec- 
tron spectrum J e (E) is well represented by 



J e {E) = 1.6xlO- 4 (£/100GeV;r 3 - 3 

(irrV^r- 1 GeV" 1 ) 



(7) 



in the energy range of 30 GeV - 1 TeV @. At a depth of 
x g cm -2 , the bremsstrahlung gamma-ray spectrum from 
the electrons with a power-law index of —3.3 is given by 



Cbrem = J e (E)C( S = 2.3) 

g— aox/Xo 



-A(s=2.3)x/X 



X 



A(s = 2.3) - a 



(8) 



where the radiation length in the atmosphere Xq = 
36.7 g cm" 2 , A(s = 2.3) = 1.674, and C(s = 2.3) = 
0.4118. The notations A and C refer to the functions 
used in electro-magnetic shower theory, defined in refer- 
ence [31|. 

The flux of the Galactic diffuse gamma-ray emission 
measured by EGRET is less than 0.1 % of the atmo- 
spheric gamma rays at depths of several g cm -2 |32| . 
Sreekumar et al. (1998) |33| reported that extra-galactic 
diffuse gamma rays observed with EGRET have a power- 
law spectrum with an index of —2.1 in the 50 MeV to 
100 GeV. Their flux at 100 GeV corresponds to - 0.5 % 
of the atmospheric gamma rays, and there are also argu- 
ments that the flux of extra-galactic diffuse gamma rays 
is lower than their derivation j34|. Therefore, we ignored 
the contribution of the flux of astrophysical gamma rays 
to obtain the atmospheric gamma-ray flux. 

The gamma-ray fluxes are normalized to 4.0 g cm 
equivalent altitude from each observation altitude of 
x g cm -2 using 



a 



a -4.0/A p 



,-4.0<t /X 



alt 



-x/Ap _ e -xa /X 



(9) 



where A p is an attenuation length of protons in the atmo- 
sphere, whose value is ~ 100 g cm~ 2 in the TeV region 
as described in section ITTT1 

The total number of observed gamma rays is 330 events 
in the energy range from 30 GeV to 8 TeV. After the 
corrections described, we derived the vertical spectrum 
of atmospheric gamma rays. Figure|S|shows the observed 
spectrum of atmospheric gamma-rays normalized at an 
altitude of 4.0 g cm~ 2 , which is well represented by 



J 7 (£) = (1.12 ±0.13) x 10- 4 (£/100GeV) 



-2.73±0.06 



(m- 2 s- 1 sr- 1 GeV^ 1 ). (10) 

The flux values and numbers of the gamma rays in each 
energy bin are listed in Table ITTT1 
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FIG. 6: The atmospheric gamma-ray spectrum observed at 
an altitude of 4.0 g cm -2 . The dash line shows the best fit 
power-law function with an index of — 2.73±0.06. 



considering the absorption of gamma-rays in the atmo- 
sphere, where <7(-E 7 ) is a production spectrum of gamma- 
rays generated in collisions between protons and atmo- 
spheric nuclei. 

In the energy range over 100 GeV, one can apply a 
scaling law in the first approximation for the production 
probability of 7r°. The ir° production rate in a single 
collision between protons and atmospheric nuclei is rep- 
resented by 



f{E„o/E p )d(E v o/E p ) 



(13) 



assuming a scaling law, where E„o is the energy of sec- 
ondary 7T° and E p the energy of parent protons. 

The production spectrum of ir° generated in collisions 
between protons with a power-law spectrum and atmo- 
spheric nuclei is now given by 



Il(E„o)dE„o = / dE p E--tf{E^/E v )d{E„o/E p ). 

JE^o 

(14) 

Putting 



TABLE III: Atmospheric gamma-ray fluxes and raw number 
of gamma rays at 4.0 g cm -2 

Energy E Raw Flux 

(GeV) (GeV) Number (m" 2 s" 1 sr" 1 GeV -1 ) 



30 
60- 
100- 
200- 
400- 
600- 
800- 
1000- 
2000- 
3000- 



50 

100 

-200 

-400 

-600 

-800 

1000 

-2000 

-3000 

-8000 



3.84X10 1 
7.68X10 1 
1.39xl0 2 
2.79xl0 2 
4.87xl0 2 
6.91xl0 2 
8.93xl0 2 
1.39 xlO 3 
2.44xl0 3 
4.76xl0 3 



8 
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(1.74±0 
(1.75±0 
(4.12±0 
(7.36±0 
(1.62±0 
(5.99±0 
(2.53±0 
(8.70±1 
(1.70±0 
(2.84±0 



)xl0~ a 
)xl0" 4 
)xl0~ 5 
)xl0~ 6 
)xl0~ 6 
)xl0~ 7 
)xl0~ 7 
)xl0~ 8 
)xl0~ 8 
)xlQ~ 9 



III. DECONVOLUTION OF PRIMARY 
PROTON SPECTRUM 

We assume that the primary protons at the top of at- 
mosphere have a power-law spectrum of 



J P (E)dE p = NE^dEp. 



(11) 



As described below, the mean free path length of 
hadronic interactions and effective attenuation length of 
protons in the atmosphere are X p ~ 80 g cm -2 and A p ~ 
100 g cm~ 2 in the TeV region, respectively [sELI^. The 
flux of atmospheric gamma rays observed at 4.0 g cm~ 2 
is represented by 



J(£ 7 ) = g(Ej 



1 exp(-4.0/A p ) - exp(-4.0cr /A" ) 



a /Xo - 1/Af 



Z p7r o = / x 1 1 f(x)dx, 
Jo 

where x = E^o / E p , we have a simple formula of 



Il(E v0 )dE„o 



Zp^aEj dE^ 



(15) 



(16) 



in which Z^o is the spectrum- weighted moment, called 
as Z factor [33 . We derived a numerical value of Z p ^o for 
the case of power-law index 7 = 2.75, using the hadronic 
interaction models Dpmjet3 ,38] and Fritiof V7.02 [39| . 
The results are shown in Fig. {7\ Here we include the 
contribution of r\ mesons to gamma rays. In implement- 
ing hadronic interaction models, we used a Monte Carlo 
simulator called COSMOS]^. Although we adopted a 
scaling-law for the deconvolution of atmospheric gamma 
rays to primary protons, there is a weak energy depen- 
dence in the Z factor of the formula (|15|) . Therefore 
we need to correct this energy dependence. As for the 
hadronic interaction models used, there are several % dif- 
ferences below 1 TeV between Dpmjet3 and Fritiof V7. 02, 
but above 1 TeV, the results from these two models agree 
well with each other. 

The production spectrum of gamma rays from a single 
7T° decay with energy E^o is given by 



dn~f 



Etj 



(17) 



Combining the formula (|16|) and (|17|) . the energy spec- 
trum of gamma rays from a single collision between pro- 
tons and atmospheric nuclei is given by 



(12) 



9(E) 



dn-y 
e dE~ t 



II(i? 7r o )dE n o 



"jiZrt-rrO 

E~ J ■ 



(18) 
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FIG. 7: Spectrum weighted moment Z p7T o of Dpmjet3 and 
Fritiof V7.02 for p + Air — ► 7T° +X with the parent proton en- 
ergy, including the contribution of r] meson decay for gamma 
rays. 



From equation (|12l) and (|18l) . the energy spectrum of 
gamma rays from interactions of primary protons in the 
atmosphere is represented by 



where 



°atm 



J 1 (E) - -Z p7r oC^ tm NE' 



1 exp(— x/A p ) — exp(— <jox/Xq) 
X p vq/Xq - 1/Ap 



(19) 



(20) 



Using this formula, we can derive the primary proton 
spectrum from the observed atmospheric gamma rays. 
The interaction length of protons in the atmosphere is 
described by the expression A p = Amp/ap 1 "' (A = 14.5), 
where A is the average mass number of atmospheric nu- 
clei, Trip is the mass of protons, and cr p " r is the cross 
section for proton-air interaction. We also derived the 
interaction length from the cross section given in ref- 
erences |2g. We derived the attenuation length 
of protons in the atmosphere using the formula A p = 
A p /(1 — Z pp — Z pn ). The spectrum- weighted moments 
of Zpp and Zp n are calculated using the hadronic inter- 
action codes of Fritiof V7.02 and Dpmjet3. We present 
the interaction length X p and attenuation length A p of 
protons for various energies in Fig. |SJ 

Here we need to include the contribution due to heavy 
primaries. We corrected for the effects of heavier nu- 
clei such as He, C, N, and O. The flux of heavier nuclei 
is ~6.4% of the flux of protons at the same energy per 
nucleon for He, and ~0.57% for C,N,Q. based on the 
JACEE and RUN JOB observations [II 111 H. Since 
interaction lengths for He and C,N,0 components are 
\ He /X p = 0.56 for He, \ N /\ P = 0.29 for C,N,0, us- 
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FIG. 8: Energy dependence of mean free path lengths of 
hadronic interactions X p and attenuation lengths A p of pro- 
tons in the atmosphere l^l |3(| . 



ing the Hagen- Watts formula ^fji atmospheric interac- 
tion contributions are C^/Cf tm = 1.73 for He, and 
Catm/Catm = 335 for C ; N >° from the formula (H1J). Ac- 
cording to Dpmjet3, the spectrum- weighted moments are 
given by Zjj^a j Z v ^a = 2.06 for He and Z^ v a / 'i? p7r o = 
4.97 for C,N,0. Therefore, using formula (HHJ, the flux of 
heavier nuclei can be converted to an equivalent proton 
flux as follows: 



^CNO =0M4 C ^ Z 



Hen 



^atm 



0.0057 



Catm Zpir 



0.32. 



(21) 



Therefore, the proton flux from the gamma rays is mul- 
tiplied by 1/1.32 = 0.76. This correction factor is almost 
same as the results of an independent-nucleon model for 
heavier nuclei as follows: 4x0.064 + 14x0.0057 = 0.34, 
that is 1/1.34 = 0.75. The uncertainty of the correction 
factor is mainly from the uncertainties of the currently 
observed flux of heavy primary 0, 0, 0] , and is esti- 
mated to be no more than ~5 %. 

We also corrected for the minor contributions of 
gamma rays from rj and K mesons, using Dpmjet3 re- 
sults. The contribution from rj mesons is ~ 0.16 as 
large as from ir° mesons. The correction for rj mesons 
is included in the calculation of Z plx a. The contribu- 
tion of gamma rays from K° mesons is ~ 0.03 as large 
as those from ir° mesons. To correct this contribution, 
we multiplied the proton flux by the correction factor of 
1/1.03 = 0.97. 

As shown from equation l|18f) . on e proton produces 
gamma-rays of ("f/2)(l/Z pK a). Hence, the energy of the 
corresponding parent protons is Ce = (2Z p7r o/7) -1 ^ 7-1 ^ 
times of that of the gamma rays. Thus we can presume 
that atmospheric gamma rays are on average produced 
from primary protons that have Ce times higher energy 
than the daughter gamma rays. Ce is ~7.2 for Dpmjet3 
with 7 = 2.75. 
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In this way, we derived the deconvolved primary proton 
spectrum as 



Jp(E p ) 



1 



a /X - 1/Ap 



7 

2 Z pw o ^ p exp(-a:/Ap) - exp(-a a:/Ao) 

x 0.76 x 0.97 x C E 7 J 7 (£ 7 ), (22) 



where £" p = CeE 7 . Using this formula, we obtained 
the primary cosmic-ray proton spectrum, using the Dp- 
mjet3 hadronic interaction model, from the atmospheric 
gamma-ray spectrum at 4.0 g cm -2 as presented in Fig. [5] 
Figure [5] shows our deconvolved proton spectrum with 
the directly observed proton spectra at high altitude ob- 
tained by other groups H 0, U d E III 111 EH Our 
proton flux with Dpmjet3 (open circles) is ~ 20 % larger 
than the extrapolation of the BESS, AMS-01, and BESS- 
TeV data. 
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FIG. 9: The proton spectra deconvolved from the atmo- 
spheric gamma-ray spectrum ; c omp ared with the spectra ob- 
served by other groups 11011 [III 111 111 El El. Open cir- 
cles show the proton spectrum derived from the gamma-ray 
spectrum with Dpmjet3. Solid circles show our final proton 
spectrum derived with 20 % larger Z^o than that of Dpmjet3. 

The uncertainties in deriving the flux of primary pro- 
tons from atmospheric gamma-ray flux come mainly from 
the hadronic interaction models, while the muon flux es- 
timated from the gamma-ray flux is less sensitive to the 
treatment of hadronic interactions, as shown in section 
II VI and appendix [BJ For hadronic interactions, there are 
uncertainties in factors such as the interaction length of 
protons in the atmosphere and the production rate of pi- 
ons. The former has an uncertainty of several % as shown 
in Fig. [HI an d the uncertain ty o f the latter is estimated to 
be much larger, ~ 20-25% |22|. Therefore, the mam sys- 
tematic error contribution in the estimation of the proton 
flux from the gamma-ray flux comes from the uncertainty 



of the Z factor for the production rate of pions. Sanuki 
et al. (2006) pointed out that the muon flux measured 
with the BESS-TeV and L3+C above 100 GeV is - 20 % 
larger than the calculated muon flux with Dpmjet3, and 
that the muon charge ratio (/i + /V - ) is inconsistent with 
that predicted by Dpmjet3. They suggested that in the 
energy range of 100 GeV to 10 TeV the production rates 
of charged pions and kaons are ~20 % larger than that 
by Dpmjet3 |2l|. They also proposed that Dpmjet3 can 
be modified in a phenomenological way, with an assump- 
tion based on the quark model, to increase by ~ 20 % 
above ~ 100 GeV the Z factors of secondary particles, 
accommodating the calculated muon charge ratio to the 
observed one. 

Referring to their modified model, we adopted 20 % 
larger Z^o than that of Dpmjet3 and derived the proton 
flux from the gamma-ray flux using formula i|22|) . The 
energy of the corresponding parent protons is Ce— 6.5 
times of that of the gamma rays. The derived proton 
spectrum (solid circles) is presented in Fig. ^\ and is well 
represented by 

J P {E) = (5.4 ± 1.2) x 10- 2 (£/100GeVr 2 - 79±0 - 06 

(m-W^GeV- 1 ) (23) 

in the energy range from 200 GeV to 50 TeV. The flux val- 
ues are also summarized in Table llVl To derive the pro- 
ton spectrum, we assumed a power-law index of —2.75. 
The fitted indexes are —2.73 for the gamma rays and 
—2.79 for the protons, slightly steeper than the gamma 
rays, because of the weak energy dependences of Z factors 
and interaction lengths, as shown in Fig. [3 and Fig. |H1 
The differences of power-law indexes from —2.75 cause 
~5 % differences for the proton flux, which are much 
smaller than the main uncertainty of hadronic interac- 
tion models of ~20 %. 

Our deconvolved proton spectrum above 10 TeV is 
somewhat smaller than the JACEE and MUBEE data, 
and agrees well with the RUNJOB data. The ATIC-2 
data show much higher flux than our proton flux. Be- 
low 10 TeV, our proton spectrum is consistent with the 
extrapolated fluxes with the BESS, AMS-01, and BESS- 
TeV data. 



IV. ESTIMATION OF HIGH-ENERGY MUON 
SPECTRUM AT HIGH ALTITUDE 

Atmospheric muons mainly come from the decay of 
charged pions and kaons. We discuss the contributions 
of pions and kaons for muon production. 

In the first approximation, the same number of ir~ and 
7r + mesons are produced as 7r°. Therefore, the production 
spectrum of charged pions can be estimated from that 
of atmospheric gamma rays through 7r°, independent of 
primary cosmic rays, such as F n ± ~ 2F^o ~ 7-F 7 , where 
F„±, F^o and F 1 are the production spectrum of 7r ± , 
7T and atmospheric gamma rays per unit depth, respec- 
tively. In more detail, the production ratio of charged 
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TABLE IV: Deconvolved primary proton fluxes 



E Flux 

(GeV) (m -2 s" 1 sr" 1 GeV -1 ) 



3.49xl0 2 


(1.52±0.38)xl0~ 


-3 


9.23xl0 2 


(9.83±2.15)xl0" 


-5 


1.84xl0 3 


(1.72±0.22)xl0~ 


-5 


3.20xl0 3 


(3.73±0.48)xl0~ 


-6 


4.52xl0 3 


(1.36±0.18)xl0~ 


-6 


5.83xl0 3 


(5.72±0.97)xl0~ 


-7 


9.04xl0 3 


(1.94±0.25)xl0~ 


-7 


1.57xl0 4 


(3.72±1.07)xl0~ 


-8 


3.03xl0 4 


(6.07±1.92)xl0~ 


-9 




Muon Energy (GeV) 



pions to neutral pions is calculated to be 1.7 using a 
hadronic interaction model of Dpmjet3. 

Since the 7r ± lifetime for decay to muons is 2.6 x 
10~ 8 sec, ir^~ mesons above ^100 GeV travel more than 
~5 km on the average and may interact with atmo- 
spheric nuclei before decay. On the other hand, charged 
kaons have larger mass, 493.6 MeV, and shorter life time, 
1.2 x 10 -8 sec. Therefore, although atmospheric muons 
are mainly produced by the decay of n^, the relative 
contribution of kaons to the muon flux increases with 
increasing energy above ^100 GeV [37| . 

Hence, on the basis of the gamma-ray flux from tt°, we 
can derive the flux of muons from ^ and K , correcting 
for muon decay and the flux of gamma rays from r\ and 
K° mesons. The muon flux is represented as the following 
formula of 

^ = ^(/.±+/*)7XT, (24) 

1 T Jrj 

where is the survival probability of muons, / ff ± is the 
ratio of the muon flux from w to the gamma-ray flux 
from 7T°, fx is a ratio of the muon flux from kaons to 
the gamma-ray flux from ir°, and is the ratio of the 
gamma-ray flux from rj + K to the gamma-ray flux from 
7T°. We derived the production ratio of K^/n^ using 
the hadronic interaction model of Dpmjet3. Figure ITUI 
shows Sp, f-K, /k, and f v as a function of muon energy. 
The detailed derivation for these parameters is given in 
appendix iBl 

We transformed the observed atmospheric gamma-ray 
spectrum to the muon spectrum using the formula (|24|l . 
and results are presented in Fig. ^] and Table [3 In 
Fig. 1111 we compare our estimated /i + + spectrum at 
4.0 g cm -2 with that of /z + and \i~ observed by BESS and 
HEAT group 0, 0] as representative values, because 
the results by the BESS, HEAT, MASS and CAPRICE 
are consistent with each other within statistical errors. 
Since the BESS and HEAT observations are performed 
at altitudes of 4.5 g cm -2 and 5.6 g cm~ 2 respectively, 
the coefficients of 4.0/4.5 and 4.0/5.6 are multiplied to 
obtain the muon flux at 4.0 g cm -2 for the compari- 
son. As shown in Fig. 1111 our flux is just on the line 



FIG. 10: Coefficients describing the relation between the 
muon flux J M and the gamma-ray flux J 1 . J M = S l _ l (f n ± + 
/k)J^/ (1 + frj), where f v = 0.16. is a survival probability 
of u . / w ±, /if, and f v are flux ratios of 7T — >/x to 7r°^7, 
K ± — >(i ± to 7T°— >7, and r],K ^'y to ty°^j, respectively. See 
text in detail. 

of extrapolation of each datum of the BESS and HEAT 
experiments and gives consistent values within statistical 
errors, although the energy region of our gamma rays is 
almost one order of magnitude higher than that of the 
directly observed muons. 

To calculate the atmospheric neutrino flux at energies 
above ~ 100 GeV with better accuracy, muon flux data 
above ~ 100 GeV measured at high altitude would be 
very useful for the calibration of the hadronic interaction 
model. However, we have no direct observations of muons 
above ~ 100 GeV at high altitude. Our estimated muon 
flux gives useful information for this purpose. 

TABLE V: Estimated atmospheric + fluxes at 
4.0 g cm" 2 



p 


Flux 


(GeV/c) 


(m" 2 s" 1 sr" 1 (GeV/c)- 1 ) 


3.85X10 1 


(1.70±0.61)xl0" a 


7.69X10 1 


(1.36±0.48)xl0~ 4 


1.39xl0 2 


(2.43±0.53)xl0~ 5 


2.79xl0 2 


(2.80±0.37)xl0~ 6 


4.87xl0 2 


(3.97±0.52)xl0~ 7 


6.91xl0 2 


(1.10±0.14)xl0~ 7 



V. SUMMARY 

We have observed atmospheric gamma rays at sev- 
eral g cm -2 altitude with emulsion chambers. Precise 
measurements of the atmospheric gamma-ray spectrum 
from hadronic interactions have been performed with the 
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FIG. 11: Comparison of muon spectra estimated by our 
gamma- ray observation with those observed by BESSflSII and 
HEATll^l. The altitudes are normalized at 4.0 g cm -2 . In 
the direct observations charged muons are measured individ- 
ually, and in our estimation all charged muons are included 
together. 



largest SfLT among all existing measurements in the en- 
ergy range 30 GeV to 8 TeV. To obtain the atmospheric 
gamma-ray flux at 4.0 g cm~ 2 from the observed data 
at each altitude, we took into account several correction 
factors such as gamma-ray detection efficiency, enhance- 
ment due to energy resolution, bremsstrahlung gamma- 
ray flux from primary electrons, and altitude conversion 
to 4.0 g cm~ 2 . Although some electronic detectors have 
uncertainty in the determination of their geometrical fac- 
tor Stt-y, that of the emulsion chambers can be estimated 
very accurately because of the precise determination of 
the track location and the simple configuration of the de- 
tector. Thus, the uncertainties of these correction factors 
and SQ-y are relatively smaller than the statistical errors. 
We deconvolved primary proton spectrum in the 200 GeV 
— 50 TeV from our gamma-ray spectrum in a reliable 
way, assuming a single interaction of each proton with 
an atmospheric nucleus. The main uncertainty in the 
deconvolution comes from hadronic interaction models, 
and we referred to a phenomenologically modified model 
in Dpmjet3 H^. While in the energy range from 100 GeV 
to 10 TeV accurate data are missing in the currently ob- 
served proton spectra, our estimated proton spectrum 
fills this gap. Our derived proton flux is consistent with 
the other observed data in the overlapping region. This 
may also indicate the validity of the hadronic interaction 
model in the TeV region proposed by Sanuki et al. |2l| . 
From the gamma-ray spectrum, we also deconvolved the 
atmospheric muon spectrum, which is consistent with di- 
rect muon observations below 10 GeV. 
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APPENDIX A: ENHANCEMENT FACTOR OF 
FLUX BY ENERGY RESOLUTION 

Defining Eq as the incident gamma-ray energy and a as 
the standard deviation of energy determination, the ob- 
served energy E has an uncertainty given by the Gaussian 
probability function 



P(E - Eo) = 



1 



■eM~(E-E ) 2 /(2a 2 )). (Al) 



2na 

As the gamma-ray spectrum is a power-law function of 
o, the enhancement factor is given by 

1 



C en h — 



E-y 



E-~<P{E - E )dE 



o- 



(A2) 



const. 



(A3) 



In the case of constant energy resolution a/Eo 
it is well represented by the series form of 

( 7 -l)( 7 -2) a 2 

Gcnh — H -Z (Tr) H > 

independent of the gamma-ray energy. The exact solu- 
tion can be presented by a hypergeometric function. In 
the case of energy resolution a/E = 15%, with 7 = 2.75, 
this enhancement factor is 1.01. In emulsion chambers, 
the energy resolution is well represented by the form of 
(|4"|. Using numerical integration with the formula I|A1|1 . 
(|A2I) . and Q, we derived the enhancement factor C cn h 
for each gamma-ray energy bin and for each emulsion 
chamber. C on h has values from 1.00 to 1.06. 



APPENDIX B: CONTRIBUTION OF PIONS AND 
KAONS FOR MUON PRODUCTION 

1. Decay factor B for each particle 

In the case of an isothermal atmosphere, secondary 
particles produced at a depth of x g cm -2 arrive at a 
depth of xo g cm -2 without decay with a survival prob- 
ability of (x/xq)^ b ^ e ^ per unit depth, neglecting the ion- 
ization loss of the particles, where B is the decay factor 
and E is the energy for each particle |4^. We define the 
decay factor B as B — Hiti/(tc), where H is the scale 
height of the atmosphere, m is particle mass, r is particle 
decay lifetime, and c is light speed. Setting scale height 
H to be 6.3 km at balloon altitude, the decay factors are 
B^ = 1.0 GeV for muons, B^ = 112.7 GeV for pions, and 
B K = 838.5 GeV for kaons. 
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2. Decay of muons 

At balloon altitudes of several g cm~ 2 residual over- 
burden, atmospheric muons are produced at a rate pro- 
portional to the transverse depth x, because the interac- 
tion length for cosmic rays on atmospheric nuclei is much 
larger than the atmospheric depth of several g cm~ 2 . 
Therefore, the survival probability of muons without 
decay is given by 



1 f x ° t 
S, = - / {-) B »' E dx 
xo Jo x o 



E 



B„ + E 



(Bl) 



3. Decay of pions 

The flux of 7T decaying per unit depth is obtained by 
multiplying B^/E into the formula ijBlf) and normalizing 
by the production rate as 



In the case of the complete decay of pions, the flux of 
muons per unit depth at a depth of a;o g cm -2 is given 
by 



E+ dE 



x 0— 5 5 

™-i - m i Je- 



■C() 



1-rv 



EI 1 



Ejr 

dE m 



xo 



1 - r f E~~< 

7 ' (B7) 



71 E n 1 - TV 0.85 7 



As a result, the ratio of the flux of muons from the 
decay of charged pions to the flux of gamma rays from 
the decay of neutral pions is given by 



1 — r~\ J E ^ B^ + E n 



(B8) 



(B2) 



Hence the flux of ir decaying until a depth of xq g cm 2 
is given by 

*o(-5^t-e;)F« = 0-®n(-5^T-Ei)Jy, ( B 3) 



' B-K + E 



V B-,1 + E 



where J 7 is the atmospheric gamma-ray flux at Xq g cm~ 2 
and given by J 7 = xqFj from the relation of F n = 
1.7 F v o — 0.857F 7 . Here, the production ratio of tt^ 1 to 
7T° is calculated to be 1.7 using the hadronic interaction 
code Dpmjet3. 

As described in reference |42( , the flux of muons pro- 
duced from the decay of pions at a depth of x g cm~ 2 , 
per unit depth, is given by 



UK 



u fi J E- * B v + E n 



where E+ and E— show the upper and lower limits of the 
energy of pions with rest mass that produce muons 
with rest mass energy E^. E+ and E— are given by 



4. Decay of kaons 

Although muons are mainly produced from the decay 
of charged pions, there is a minor contribution from kaons 
[3?| . While there are K^, Kg and -ftf£, only the charged 
kaons decay directly to muons via K ± — > fx^ + v. For 
the production of muons the contribution of kaons can 
be treated in the same way as pions. We derive the rela- 
tion between kaons and gamma rays through the produc- 
tion ratio of if ± /7r ± . The production ratio of K^/tt^ is 
0.14 in the parent proton energy of 200 GeV — 40 TeV 
according to the hadronic interaction code in Dpmjet3. 
The decay mode K ± ^fi ± dominates (63 % of decays) 
among the charged kaon decay mode [43j. As a result, 
we can obtain the ratio of the flux of muons from kaons 
to the flux of gamma rays from neutral pions as 



f K (E, t ) = 0.14x0.63x 



0.85 7 ^ 



1 

E^/tk 



K 

B 



K 



B 



K 



E 



K 



-)E 



,dE 



K 



K 



E 



K 



(B9) 



and 



* „2 



E+ = ( — )( 



w.fE^+p^c 



{'^?E, 



E„ 



(B5) where tk is (m^/rriK) 2 — 0.0459. Although there is a 
decay mode (21 % of decays), the contribution 



£ =( ^y )2 , (B6) 
m M m M c 2 



where r n is (m M / 'm^) 2 — 0.5733, is momentum of 
pion, and the notation V means a static system. 



± is less than ~1 % for the total muon 



if* 

Of K ± — > 7T^ — » fl 

flux. 

Since kaons have larger rest mass and shorter life time, 
the decay factor of kaons, Bk = 838.5 GeV, becomes 
larger than that of pions, B^ = 112.7 GeV. Therefore, 
the contribution of kaons continues into a higher energy 
region than pions. 
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